1 Moran’s I vs window size

stromaLabels = c("CD20pos Bcells", "CD4posTcell", "CD8posTcell", "CollagenCD45", "Stromal something", "Treg", "Macrophage")
epiLabels = levels(sc_HTAfull$label)[ !levels(sc_HTAfull$label) %in% stromaLabels]
sc_HTAfull$epiCell = ifelse(sc_HTAfull$label %in% stromaLabels, 0, 1)
sc_HTAfull[, levels(sc_HTAfull$label)] = model.matrix( ~ -1 + label, data=sc_HTAfull)
dat <- filter(sc_HTAfull, epiCell==1)
dat <- filter(dat, region=="region_001")
dat <- filter(dat, region=="region_001")
slides <- unique(dat$SlideID)[1:10]
dat <- filter(dat, SlideID %in% slides)
ADSSL <- c(2,3,5,7,9,1,4,6,8,10)
rs <- c(50, 100, 250, 500, 600)
markers <- c("TumorEpiMuc2","TumorStem","IEL","Goblet")

## plot function: Moran's I
moransI.df.transfer <- function(result.df){
  result.df$marker1 <- rep(c("TumorEpiMuc2","TumorStem","IEL","Goblet"), each=4)
  result.df$marker2 <- rep(c("TumorEpiMuc2","TumorStem","IEL","Goblet"), 4)
  plot.df <- melt(result.df,variable.name = "radii", id.vars = c("marker1","marker2","slide","TumorType"), 
                  value.name = "MoransI")
  plot.df$radii <- as.numeric(levels(plot.df$radii))[plot.df$radii]
  plot.df
}
plot.moransI <- function(result.df){
  plot.df <- moransI.df.transfer(result.df)
  p <- ggplot(plot.df, aes(x=radii, y=MoransI, group=slide)) +
  geom_line(aes(color=TumorType), alpha=0.5) +
  geom_point(shape=21, color="black", size=0.5, alpha=0.3)+ theme_bw()+
  facet_grid(rows=vars(marker1), cols=vars(marker2), scales = "free_y")
  p
}
# plot function for noramlized moran's I
mI.standard.transfer <- function(result.df){
  result.df$marker <- c("TumorEpiMuc2","TumorStem","IEL","Goblet")
  plot.df <- melt(result.df,variable.name = "radii", id.vars = c("marker","slide","TumorType"), value.name = "MoransI")
  plot.df$radii <- as.numeric(levels(plot.df$radii))[plot.df$radii]
  plot.df
}
plot.moransI.standard <- function(result.df){
  plot.df <- mI.standard.transfer(result.df)
  p <- ggplot(plot.df, aes(x=radii, y=MoransI, group=slide)) +
  geom_line(aes(color=TumorType), alpha=0.5) +
  geom_point(shape=21, color="black", size=0.5, alpha=0.3)+ theme_bw()+
  facet_wrap(~marker)
  p
}
# plots_1 <- plots_2 <- list()
result1.all <- result2.all <-result3.all <-result4.all <-list()
for(j in ADSSL){
  result1.df <- data.frame(matrix(NA, ncol=length(rs), nrow=16)) # moran's I
  result2.df <- data.frame(matrix(NA, ncol=length(rs), nrow=4)) # Moran's I standardized
  colnames(result1.df) <- colnames(result2.df) <- as.character(rs)
  result1.df$slide <- result2.df$slide <-  slides[j]
 
  
  result3.list <- result4.list <- list()
  for(i in 1:length(rs)){
    result3.df <- data.frame(matrix(NA, ncol=4, nrow=sum(dat$SlideID==slides[j]))) # LISA
    result4.df <- data.frame(matrix(NA, ncol=4, nrow=sum(dat$SlideID==slides[j]))) # LISA standardized
    colnames(result3.df) <- colnames(result4.df) <-markers
    result3.df$slide <- result4.df$slide <- slides[j]
    
    dat_i <- filter(dat, SlideID == slides[j])
    wlist = multiplexMoran::coords2listw(dat_i[,c("x", "y")], k = 500, maxDist = rs[i] )
    moranMat = multiplexMoran::moran(select(dat_i, TumorEpiMuc2,TumorStem,IEL,Goblet),
    listw = wlist, demean=TRUE)
    result1.df[,i] <- as.vector(moranMat)
    for (k in 1:length(markers)){
       z_hat <- (moran.test(dat_i[, markers[k]], listw=wlist)$estimate)[1]
       result2.df[k,i] <- z_hat
       
       lisaVec = spdep::localmoran(dat_i[, markers[k]], listw = wlist)
       result3.df[,k] <- lisaVec[,1]
       result4.df[,k] <- lisaVec[,4]
    }
    result3.df$r <- result4.df$r <- rs[i]
    result3.list[[i]] <- result3.df
    result4.list[[i]] <- result4.df
  }
result1.all[[j]] <- result1.df
result2.all[[j]] <- result2.df
result3.all[[j]] <- result3.list
result4.all[[j]] <- result4.list
  # plots_1[[j]] <- plot.moransI(result1.df, name.slide = slides[j])
  # plots_2[[j]] <- plot.moransI.standard(result2.df, name.slide = slides[j])
}
moran_result_list <- list(result1.all, result2.all, result3.all, result4.all)
save(moran_result_list, file="moran_0718.RData")
setwd("/home/xiongj3/hist_fit_moran/atlasAnalysis")
load("moran_0718.RData")
# for( i in 1:5 ){
#   print(plots_1[[ADSSL[i]]])
# }
result1.all.df <- do.call("rbind", moran_result_list[[1]])
result1.all.df$TumorType <- rep(c("Adenoma", "SSL"), each=(nrow(result1.all.df)/2))
plot.moransI(result1.all.df)

# for( i in 6:10){
#   print(plots_1[[ADSSL[i]]])
# }

2 Standardized Moran’s I vs window size

# for( i in 1:5 ){
#   print(plots_2[[ADSSL[i]]])
# }
result2.all.df <- do.call("rbind", moran_result_list[[2]])
result2.all.df$TumorType <- rep(c("Adenoma", "SSL"), each=(nrow(result2.all.df)/2))
plot.moransI.standard(result2.all.df)

# for( i in 6:10 ){
#   print(plots_2[[ADSSL[i]]])
# }

3 Number of cells v.s. Index

Effect of number of cells on value of Moran’s I

library(dplyr)
cell.num <- as.data.frame(table(dat$SlideID))
colnames(cell.num) <- c("slide", "cell.count")
plot.df <- moransI.df.transfer(result1.all.df)
cvi.plot.df <- left_join(plot.df, cell.num, by="slide")
ggplot(cvi.plot.df, aes(x=cell.count, y=MoransI))+
  geom_point(aes(color=TumorType), alpha=0.5)+facet_wrap(~radii, nrow=1)+theme_bw()

Effect of number of cells on value of standardized Moran’s I

plot.df <- mI.standard.transfer(result2.all.df)
cvi.plot.df <- left_join(plot.df, cell.num, by="slide")
ggplot(cvi.plot.df, aes(x=cell.count, y=MoransI))+
  geom_point(aes(color=TumorType), alpha=0.5)+facet_wrap(~radii, nrow=1)+theme_bw()

4 LISA

LISA is taken \(\sqrt[\leftroot{-1}\uproot{2}\scriptstyle 5]{LISA}\qquad\) for illustration purposes

result3.list <- moran_result_list[[3]]
adssl <- rep(c("Adenoma", "SSL"), each=5)
library(cowplot)
cube.root <- function(number){
  sign <- (-1)^(1*(number<0))
  sign*abs(number)^(1/5)
}
plot.lisa <- function(plot.lisaMat, root.trans=TRUE){
  plot.lisa.df <- melt(plot.lisaMat,variable.name = "Markers", id.vars = c("slide","r","x","y"), 
                  value.name = "LISA")
  if(root.trans==TRUE){
    plot.lisa.df$LISA <- cube.root(plot.lisa.df$LISA)
  }
  
  ggplot(data=plot.lisa.df, aes(x=x, y=y, color=LISA)) +
    geom_point(size=0.05, alpha=0.3)+
    scale_color_gradient2()+
    theme_dark()+theme(legend.position="top",legend.text = element_text(angle = 50))+
    facet_grid(rows=vars(r), cols=vars(Markers), scales = "free_y")
}
lisa.df.print <- function(i){
  result3.df <- result3.list[[i]]
  plot.lisaMat <- do.call("rbind", result3.df)
  plot.lisaMat$x <- dat[dat$SlideID==plot.lisaMat$slide[1],"x"]
  plot.lisaMat$y <- dat[dat$SlideID==plot.lisaMat$slide[1],"y"]
  p1 <- plot.lisa(plot.lisaMat)+ggtitle(paste(slides[i], adssl[i]))
  dat.i <- dat[dat$SlideID==plot.lisaMat$slide[1],]
  dat.i <- subset(dat.i, dat.i$label %in% markers)
  dat.i$label <- factor(dat.i$label, levels = markers)
  p2 <- ggplot(dat.i, aes(x=x, y=y))+geom_point(size=0.05)+
    theme_bw()+facet_wrap(~label, scales = "fixed", nrow=1)
  pp <- plot_grid(p1, p2, nrow = 2, rel_heights = c(4/5, 1/5))
  print(pp)
}

4.1 HTA11_10167_0000_01_01

lisa.df.print(1)

4.2 HTA11_10623_0000_01_01

lisa.df.print(2)

4.3 HTA11_10711_0000_01_01

lisa.df.print(3)

4.4 HTA11_4255_0000_02_02

lisa.df.print(4)

4.5 HTA11_6298_0000_04A_03

lisa.df.print(5)

4.6 HTA11_6801_0000_01_01

lisa.df.print(6)

4.7 HTA11_7179_0000_02_02

lisa.df.print(7)

4.8 HTA11_7663_0000_01_01

lisa.df.print(8)

4.9 HTA11_7862_0000_02_02

lisa.df.print(9)

4.10 HTA11_7956_0000_02_05

lisa.df.print(10)

5 Standardized Lisa

result4.list <- moran_result_list[[4]]
z.lisa.df.print <- function(i){
  result4.df <- result4.list[[i]]
  plot.lisaMat <- do.call("rbind", result4.df)
  plot.lisaMat$x <- dat[dat$SlideID==plot.lisaMat$slide[1],"x"]
  plot.lisaMat$y <- dat[dat$SlideID==plot.lisaMat$slide[1],"y"]
  p1 <- plot.lisa(plot.lisaMat, root.trans=FALSE)+
    ggtitle(paste(slides[i], adssl[i]))
  dat.i <- dat[dat$SlideID==plot.lisaMat$slide[1],]
  dat.i <- subset(dat.i, dat.i$label %in% markers)
  dat.i$label <- factor(dat.i$label, levels = markers)
  p2 <- ggplot(dat.i, aes(x=x, y=y))+geom_point(size=0.05)+
    theme_bw()+facet_wrap(~label, scales = "fixed", nrow=1)
  pp <- plot_grid(p1, p2, nrow = 2, rel_heights = c(4/5, 1/5))
  print(pp)
}

5.1 HTA11_10167_0000_01_01

z.lisa.df.print(1)

5.2 HTA11_10623_0000_01_01

z.lisa.df.print(2)

5.3 HTA11_10711_0000_01_01

z.lisa.df.print(3)

5.4 HTA11_4255_0000_02_02

z.lisa.df.print(4)

5.5 HTA11_6298_0000_04A_03

z.lisa.df.print(5)

5.6 HTA11_6801_0000_01_01

z.lisa.df.print(6)

5.7 HTA11_7179_0000_02_02

z.lisa.df.print(7)

5.8 HTA11_7663_0000_01_01

z.lisa.df.print(8)

5.9 HTA11_7862_0000_02_02

z.lisa.df.print(9)

5.10 HTA11_7956_0000_02_05

z.lisa.df.print(10)